home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / WINDOWS / PROFFT.ARJ / FFT.CPP < prev    next >
C/C++ Source or Header  |  1992-04-27  |  4KB  |  150 lines
  1. /****************************************************************************
  2.     FFT.CPP            Denne filen inneholder (i denne rekkef°lgen):
  3. *****************************************************************************
  4. #define M_PI 3.14159265358979323846
  5. int sincos_length = 0;
  6. int newpow = 0;
  7. prcomplex *sincos_tab;
  8. int sincos(int length)
  9. int fft(prcomplex *z1, prcomplex *z2, int n)
  10. ****************************************************************************/
  11. #include <math.h>
  12. #include <stdio.h>
  13. #include <stdlib.h>
  14.  
  15. #pragma hdrstop
  16. #ifndef __PRCOMP_H
  17. #include "prcomp.h"
  18. #endif
  19.  
  20. //  Definerer M_PI som brukes ved oppretting av sinus cosinus tabellen
  21. #define M_PI 3.14159265358979323846
  22.  
  23. //  Husker lengden pσ tabellen
  24. int sincos_length = 0;
  25. int newpow = 0;
  26.  
  27. //  Peker til tabellen
  28. prcomplex *sincos_tab;
  29.  
  30. int sincos(int length)
  31. /****************************************************************************
  32. Denne lager en tabell for raskt oppslag pσ komplekse sinus og cosinusverdier
  33. som brukes av fft algortimen. Denne er sakset fra et dokument som ble delt
  34. ut i forbindelse med oppdraget og stammer fra Universitetet i Oslo.
  35.  
  36. int length        Lengden pσ tabellen.
  37.  
  38. Returnerer:        En toerpotens av st°rrelsen (?!?).
  39.  
  40. Kodet av:   MK    02.04.92  Modifisert for vσre komplekse talltyper.
  41. *****************************************************************************/
  42. {
  43.     int i, totlen;
  44.  
  45.     //  Hvis tabellen er opprettet allerede
  46.     if (sincos_length == length) return(newpow);
  47.  
  48.     //  Hvis lengden er ulik 0 men ikke lik lengden sσ fjern tabellen
  49.     //  fra minnet.
  50.     if (sincos_length != 0) free(sincos_tab);
  51.  
  52.     //  Sett riktige startverdier
  53.     sincos_length = 0;
  54.     i = 1; newpow = 0;
  55.  
  56.     while (i < length)
  57.     {
  58.         i <<= 1;
  59.         newpow++;
  60.     }
  61.     if (i != length)
  62.         newpow = 0;
  63.     totlen = length << 1;
  64.  
  65.     //  Alloker minne til tabellen
  66.     sincos_tab = (prcomplex*) malloc(sizeof(prcomplex)*(totlen+1));
  67.     if (sincos_tab==NULL)
  68.         return(-1);
  69.  
  70.     //  Ta vare pσ lengden av tabellen til senere
  71.     sincos_length = length;
  72.     for (i=0; i<=totlen; i++)
  73.     {
  74.         //  Fyll inn de kalkulerte verdiene i tabellen.
  75.         sincos_tab[i].re = (real)cos((double) i / (double) length * M_PI);
  76.         sincos_tab[i].im = (real)-sin((double) i / (double) length * M_PI);
  77.     }
  78.     return(newpow);
  79. }
  80.  
  81. int fft(prcomplex *z1, prcomplex *z2, int n)
  82. /****************************************************************************
  83. Dette er en rask FFT algoritme sakset fra et dokument som ble delt
  84. ut i forbindelse med oppdraget og stammer fra Universitetet i Oslo. Den
  85. er derfor ikke dokumentert i noen sµrlig grad.
  86.  
  87. prcomplex *z1        Peker til f°rste arbeidsbuffer. Mσ vµre ferdig allokert
  88.                                 og utfyllt ved kall.
  89.  
  90. prcomplex *z2   Peker til andre arbeidsbuffer. Mσ vµre ferdig allokert.
  91.  
  92. Returnerer:        1  Hvis f°rste arbeidsbuffer inneholder transformen.
  93.                             2  Hvis andre arbeidsbuffer inneholder transformen.
  94.                             0  Hvis transformen feiler.
  95.  
  96. Kodet av:   MK    02.04.92  Modifisert for vσre komplekse talltyper.
  97.                         MK  03.04.92  Fjernet un°dvendig kode.
  98. *****************************************************************************/
  99. {
  100.     int power, inda, nh, dir;
  101.     int zstep, inzee, ind, ia;
  102.     prcomplex *ar1, *ar2;
  103.     prcomplex arg2;
  104.     prcomplex *zind, *fpt, *tpt, *pt;
  105.     real fn;
  106.  
  107.     //  Hvis tabellen ikke kunne lages.
  108.     if (sincos(n)<=0) return(0);
  109.  
  110.     nh = n / 2;
  111.     power = dir = 1;
  112.     ar1 = z1;
  113.     ar2 = z2;
  114.  
  115.     for (zstep=n;zstep!=1;zstep>>=1)
  116.     {
  117.         zind = &sincos_tab[n];
  118.         ia = power;
  119.         while (ia--)
  120.         {
  121.             zind = &zind[-zstep];
  122.             fpt = &ar1[inda=nh+ia-power];
  123.             tpt = &ar2[2*inda-ia];
  124.             while (inda >= 0)
  125.             {
  126.                 pt = &fpt[nh];
  127.                 arg2.re = zind->re * pt->re - zind->im * pt->im;
  128.                 arg2.im = zind->re * pt->im + zind->im * pt->re;
  129.                 pt = &tpt[power];
  130.                 tpt->re = fpt->re + arg2.re;
  131.                 tpt->im = fpt->im + arg2.im;
  132.                 pt->re = fpt->re - arg2.re;
  133.                 pt->im = fpt->im - arg2.im;
  134.  
  135.                 fpt = &fpt[-power];
  136.                 tpt = &tpt[-power];
  137.                 tpt = &tpt[-power];
  138.                 inda -= power;
  139.             }
  140.         }
  141.         pt = ar1;
  142.         ar1 = ar2;
  143.         ar2 = pt;
  144.         dir = 3 - dir;
  145.         power <<= 1;
  146.     }
  147.     return(dir);
  148. }
  149.  
  150.